function diff=find_q_with_inc(p_interest,q_interest,y_interest,tau,theta,config)

K1=config.K1;
K2=config.K2;
K3=config.K3;
p_lb=config.p_lb;
p_ub=config.p_ub;
q_lb=config.q_lb;
q_ub=config.q_ub;
y_lb=config.y_lb;
y_ub=config.y_ub;

mu1=0;
sigma1=0*p_interest; 
mu2=NaN;
mu3=NaN;
sigma2=NaN;
sigma3=NaN;
fdist1='Normal';
fdist2='';
fdist3='';
alpha1=NaN;
alpha2=NaN;
alpha3=NaN;

[TT_interest, ~, ~, ~] = construct_chebyshev_berkson_Ponly_twodim_inc(p_interest,q_interest,y_interest,K1,K2,K3,p_lb,q_lb,y_lb,p_ub,q_ub,y_ub,NaN,NaN,NaN,NaN,NaN,NaN,mu1,mu2,mu3,sigma1,sigma2,sigma3,fdist1,fdist2,fdist3,alpha1,alpha2,alpha3);

Ginv_hat=TT_interest*theta;
diff=tau-Ginv_hat;

	%





